pcvrDDPSC Data Science Core, August 2023
Josh Sumner
pcvr GoalsgrowthSSfitGrowthbrms directlypcvr GoalsCurrently pcvr aims to:
There is room for goals to evolve based on feedback and scientific needs.
Pre-work was to install R, Rstudio, and pcvr with dependencies.
plantCV allows for user friendly high throughput image based phenotyping
Resulting data follows individuals over time, which changes our statistical needs.
Longitudinal Data is:
Bayesian modeling allows us to account for all these problems via a more flexible interface than frequentist methods.
Bayesian modeling also allows for non-linear, probability driven hypothesis testing.
In a Bayesian context we flip “random” and “fixed” elements.
| Fixed | Random | Interpretation | |
|---|---|---|---|
| Frequentist | True Effect | Data | If the True Effect is 0 then there is an \(\alpha\cdot100\)% chance of estimating an effect of this size or more. |
| Bayesian | Data | True Effect | Given the estimated effect from our data there is a P probability of the True Effect being a difference of at least X |
There are 6 main growth models supported in pcvr.
3 are asymptotic, 3 are non-asympototic.
There are also two double sigmoid curves intended for use with recovery experiments.
growthSSAny of the models shown above can be specified easily using growthSS.
growthSS - formWith a model specified a rough formula is required to parse your data to fit the model.
The layout of that formula is:
outcome ~ time|individual/group
growthSS - form 2Here we would use y~time|id/group
growthSS - form 3Generally it makes sense to visually check that your formula covers your experimental design.
Note that it is fine for id to be duplicated between groups, but not within groups
growthSS - sigmaRecall the heteroskedasticity problem, shown again with our simulated data:
growthSS - sigmaThere are lots of ways to model a trend like that we see for sigma.
pcvr offers four options through growthSS.
growthSS - Intercept sigma| Pros | Cons |
|---|---|
| Faster model fitting | Very inaccurate intervals at early timepoints |
growthSS - Linear sigma| Pros | Cons |
|---|---|
| Models still fit quickly | Variance tends to increase non-linearly |
| Easy testing on variance model |
growthSS - Linear sigmagrowthSS - Gompertz sigma| Pros | Cons |
|---|---|
| Models fit much faster than splines | Slightly slower than linear sub-models |
| Variance is often asymptotic | Requires priors on sigma model |
| Easy testing on variance model |
growthSS - Gompertz sigmagrowthSS - Spline sigma| Pros | Cons |
|---|---|
| Very flexible and accurate model for sigma | Significantly slower than other options |
| Fewer priors | Splines can be a black-box |
growthSS - Spline sigmagrowthSS - other sigma modelsYou can always add a new sigma formula if something else fits your needs better.
growthSS - priorsBayesian statistics combine prior distributions and collected data to form a posterior distribution.
Luckily, in the growth model context it is pretty easy to set “good priors”.
growthSS - priors“Good priors” are generally mildly informative, but not very strong.
They provide some well vetted evidence, but do not overpower the data.
growthSS - priorsFor our setting we know growth is positive and we should have basic impressions of what sizes are possible.
At the “weakest” side of these priors we at least know growth is positive and the camera only can measure some finite space.
growthSS - priors 2Default priors in growthSS are log-normal
\(\text{log}~N(\mu, 0.25)\)
This has the benefit of giving a long right tail and strictly positive values while only requiring us to provide \(\mu\).
growthSS - priors 3We can see what those log-normal distributions look like with plotPrior.
growthSS - priors 3growthSS - priors 4Those distributions can still be somewhat abstract, so we can simulate draws from the priors and see what those values yield in our growth model.
growthSS - priors 4growthSS - priors 5Our final call to growthSS will look like this for our sample data.
fitGrowthNow that we have the components for our model from growthSS we can fit the model with fitGrowth.
fitGrowth 2This will call Stan outside of R to run Markov Chain Monte Carlo (MCMC) to get draws from the posterior distributions. We can control how Stan runs with additional arguments to fitGrowth, although the only required argument is the output from growthSS.
fitGrowth 2Here we specify our ss argument to be the output from growthSS and tell the model to use 4 cores so that the chains run entirely in parallel, but the rest of this model is using defaults.
fitGrowth 2Note that there are lots of arguments that can be passed to brms::brm via fitGrowth.
One that can be very helpful for fitting complex models is the control argument, where we can control the sampler’s behavior.
fitGrowth 2adapt_delta and tree_depth are both used to reduce the number of “divergent transitions” which are times that the sampler has some departure from the True path and which can compromise the results.
fitGrowth 3fitGrowth returns a brmsfit object, see ?brmsfit and methods(class="brmsfit") for general information.
Within pcvr there are several functions for visualizing these objects.
brmPlot can be used to plot credible intervals of your model.
These plots can show one of the benefits of an asymptotic sub model well.
Here we check our model predictions to 35 days.
And now we check those predictions from a spline model, where the basis functions were not expecting to have data past day 25.
We can also plot the posterior distributions and test hypotheses with brmViolin.
Here hypotheses are tested with brms::hypothesis.
brms::hypothesis allows for incredibly flexible hypothesis testing.
Here we test for an asymptote for group A at least 20% larger than that of group B.
Here we test the model of our bellwether data to see if the ratio of two genotypes is 10% larger when fertilizer is added.
Hypothesis Estimate
1 ((A_group0.B73/A_group0.MM))-(1.1*(A_group50.B73/A_group50.MM)) > 0 0.2412093
Est.Error CI.Lower CI.Upper Evid.Ratio Post.Prob Star
1 0.3224859 -0.2874997 0.7497439 3.728132 0.7885
brms directlyThese functions are all to help use common growth models more easily.
The choices in pcvr are a tiny subset of what is possible with brms, which itself is more limited than Stan.
brms directlyOur gompertz sigma model looks like this in brms:
prior1 <- prior(gamma(2,0.1), class="nu", lb=0.001)+
prior(lognormal(log(130), .25),nlpar = "A", lb = 0) +
prior(lognormal(log(12), .25), nlpar = "B", lb = 0) +
prior(lognormal(log(1.2), .25), nlpar = "C", lb = 0)+
prior(lognormal(log(25), .25),nlpar = "subA", lb = 0) +
prior(lognormal(log(20), .25), nlpar = "subB", lb = 0) +
prior(lognormal(log(1.2), .25), nlpar = "subC", lb = 0)
form_b <- bf(y ~ A*exp(-B*exp(-C*time)),
nlf(sigma ~ subA*exp(-subB*exp(-subC*time))),
A+B+C+subA+subB+subC ~ 0+group,
autocor = ~arma(~time|sample:group,1,1),
nl = TRUE )
fit_g2 <- brm(form_b, family = student, prior = prior1, data = simdf,
iter = 1000, cores = 4, chains = 4, backend = "cmdstanr", silent = 0,
control = list(adapt_delta = 0.999,max_treedepth = 20),
init = 0 ) # chain initialization at 0 for simplicitybrms directlyIt can be more work to try new options in brms or Stan, but if you have a situation not well represented by the existing models then it may be necessary.
If you run into a novel situation please reach out and we will try to come up with a solution and add it to pcvr if possible.
Good ways to reach out are the help-datascience slack channel and pcvr github repository.